%main replication file for JPE Macro article "Price Setting with Customer Capital: Sales, Teasers and Rigidity" by Leena Rudanko
%this file produces figures 1 and 2 
clear all;close all;

options = optimoptions('fsolve','TolFun',1e-9,'TolX',1e-9,'StepTolerance',1e-9,'Display','final');

%parameters
be=.95^(1/12);d=1-(1-.07)^(1/12);del=1-(1-.4)^(1/12);u=1;g=1.022^(1/12)-1;ent=1-(1+g)*(1-d);nu=ent/d;kap0=1;mu0=1;
disc_target=.23;markup_target=1.15;al_target=.22;%note: resulting values for model with shocks provided below
[eta0,phi,c] = cal_fun_find_pscc(disc_target,markup_target,al_target,d,del,g);

%discretize AR(1) with rhoz,sigz based on Rouwenhorst method
n=5;rhoz=0;pz=(1+rhoz)/2;
M2=[pz,1-pz;1-pz,pz];
M3=pz*[M2,zeros(2,1);zeros(1,2),0]+(1-pz)*[zeros(2,1),M2;0,zeros(1,2)]+(1-pz)*[zeros(1,2),0;M2,zeros(2,1)]+pz*[0,zeros(1,2);zeros(2,1),M2];
M3=diag([1,.5,1])*M3;
M4=pz*[M3,zeros(3,1);zeros(1,3),0]+(1-pz)*[zeros(3,1),M3;0,zeros(1,3)]+(1-pz)*[zeros(1,3),0;M3,zeros(3,1)]+pz*[0,zeros(1,3);zeros(3,1),M3];
M4=diag([1,.5,.5,1])*M4;
M5=pz*[M4,zeros(4,1);zeros(1,4),0]+(1-pz)*[zeros(4,1),M4;0,zeros(1,4)]+(1-pz)*[zeros(1,4),0;M4,zeros(4,1)]+pz*[0,zeros(1,4);zeros(4,1),M4];
M5=diag([1,.5,.5,.5,1])*M5;
Pzz=M5;%transition matrix

%produce figure 1: temporary sale pricing

%find steady state
shat_fun=@(t) (eta_el_pscc(t,mu0,eta0)/(1-eta_el_pscc(t,mu0,eta0))*phi/(phi-1)-1)/eta_pscc(t,mu0,eta0);
al_fun=@(t) ((1-ent)/((1-d)*(1-del))-1)/(eta_pscc(t,mu0,eta0)*shat_fun(t));

res_fun=@(t) (u-c-be*(1-ent-(1-d)*(1-del))*kap_pscc(shat_fun(t),kap0,phi)/(eta_pscc(t,mu0,eta0)*shat_fun(t))-(1-be+be*ent)*kap_s_pscc(shat_fun(t),kap0,phi)/eta_pscc(t,mu0,eta0))/(-kap_n_pscc(shat_fun(t),kap0,phi)*(1-be*(1-d)*(1-del)*(1-mu_pscc(t,mu0,eta0))))-1;
%note: need to look for eqm in the range where the eqns determine it(see proof)
res_fun_num=@(t) (u-c-be*(1-ent-(1-d)*(1-del))*kap_pscc(shat_fun(t),kap0,phi)/(eta_pscc(t,mu0,eta0)*shat_fun(t))-(1-be+be*ent)*kap_s_pscc(shat_fun(t),kap0,phi)/eta_pscc(t,mu0,eta0));
th_lowerbound=fsolve(res_fun_num,1,options);
th_upperbound=(phi/(phi-1))^(1/eta0);
th_ss=fsolve(res_fun,(th_lowerbound+th_upperbound)/2,options);
shat_ss=shat_fun(th_ss);

%set value of lambda, confirm below that coincides with stationary equilibrium
la=0.0186;

%equilibrium takes form of one of several cases, which differ in which productivity states firms actively sell (or randomize) 
%consider six cases explicitly (calibration implies that firms actively selling in only limited number of states)

%initial guesses for each case use steady state values
shat_guess6=[shat_ss;shat_ss;shat_ss];th_guess6=[th_ss;th_ss;th_ss];
shat_guess5=[shat_ss;shat_ss];th_guess5=[th_ss;th_ss];al_guess5=.5;
shat_guess4=[shat_ss;shat_ss];th_guess4=[th_ss;th_ss];
shat_guess3=shat_ss;th_guess3=th_ss;al_guess3=.5;
shat_guess2=shat_ss;th_guess2=th_ss;
al_guess1=.5;
th_guess_rand=th_ss;
M=inv(eye(n)-be*(1-d)*(1-del)*Pzz);

%can find values for randomization state separately (independent of productivity)
shat_fun=@(t) (eta_el_pscc(t,mu0,eta0)/(1-eta_el_pscc(t,mu0,eta0))*phi/(phi-1)-1)/eta_pscc(t,mu0,eta0);
res_fun = @(t) la-mu_pscc(t,mu0,eta0)*(-kap_n_pscc(shat_fun(t),kap0,phi));
[th_rand,fval,exit_flag_rand]= fsolve(res_fun,th_guess_rand,options);
shat_rand=shat_fun(th_rand);
S_rand=kap_s_pscc(shat_rand,kap0,phi)/eta_pscc(th_rand,mu0,eta0)+la/mu_pscc(th_rand,mu0,eta0);

%continuation in standard deviation of shock
vol_vec=.001:.001:.10;
case_ind=0;
for i=1:length(vol_vec),i/length(vol_vec)
    %productivity vector
    sigz=vol_vec(i);
    phiz=sigz*sqrt(n-1);
    z=ones(n,1)+[phiz;phiz/2;0;-phiz/2;-phiz];
    cz=c*flip(z);%cost vector

    if case_ind==0||case_ind==6 %case 6: all firms in top three productivity states actively selling, no one else actively selling
        
        %find shat's satisfying equilibrium conditions, given lambda 
        res_fun = @(x) resid_dim5case6_pscc(x,la,u,mu0,phi,kap0,eta0,be,del,d,Pzz,cz,th_guess6,M);
        [x_6,fval,exit_flag] = fsolve(res_fun,shat_guess6,options);exit_flag 
        shat_1=x_6(1);shat_2=x_6(2);shat_3=x_6(3);
        
        %find corresponding theta's
        th_res = @(t) la-mu_pscc(t,mu0,eta0)*eta_el_pscc(t,mu0,eta0)/(1-eta_el_pscc(t,mu0,eta0))*shat_1*kap_s_pscc(shat_1,kap0,phi)/(1+eta_pscc(t,mu0,eta0)*shat_1);
        [th_1,fval,exit_flag1] = fsolve(th_res,th_guess6(1),options);
        th_res = @(t) la-mu_pscc(t,mu0,eta0)*eta_el_pscc(t,mu0,eta0)/(1-eta_el_pscc(t,mu0,eta0))*shat_2*kap_s_pscc(shat_2,kap0,phi)/(1+eta_pscc(t,mu0,eta0)*shat_2);
        [th_2,fval,exit_flag2] = fsolve(th_res,th_guess6(2),options);
        th_res = @(t) la-mu_pscc(t,mu0,eta0)*eta_el_pscc(t,mu0,eta0)/(1-eta_el_pscc(t,mu0,eta0))*shat_3*kap_s_pscc(shat_3,kap0,phi)/(1+eta_pscc(t,mu0,eta0)*shat_3);
        [th_3,fval,exit_flag3]= fsolve(th_res,th_guess6(3),options);
        
        %update initial guesses
        shat_guess6=[shat_1;shat_2;shat_3];th_guess6=[th_1;th_2;th_3];
        
        %firm value from active selling versus not in state three
        test_3=-kap_n_pscc(shat_3,kap0,phi)-la/mu_pscc(th_3,mu0,eta0);
        
        if test_3>=0, case_ind=6; disp('case 6');%if prefer active selling in state three, label as case 6 (calibration implies active selling relatively scarce) 
            %produce corresponding transition matrix and value vectors
            A=[eta_pscc(th_1,mu0,eta0)*shat_1+1,0,0,0,0;0,eta_pscc(th_2,mu0,eta0)*shat_2+1,0,0,0;0,0,eta_pscc(th_3,mu0,eta0)*shat_3+1,0,0;0,0,0,1,0;0,0,0,0,1];
            shat=[shat_1;shat_2;shat_3;0;0];th=[th_1;th_2;th_3;0;0];
            al=[1;1;1;0;0];%probabilities of active selling across states
            Pr=(u-be*(1-d)*(1-del)*la)/(1-be*(1-d)*(1-del))*ones(n,1);
            for ii=1:n, Ps(ii,1)=Pr(ii,1)-la/mu_pscc(th(ii,1),mu0,eta0);end
            pr=Pr-be*(1-d)*(1-del)*Pzz*(Pr+diag(al)*(Ps-Pr));ps=Ps-be*(1-d)*(1-del)*Pzz*(Pr+diag(al)*(Ps-Pr));
        else case_ind=5;%otherwise move to consider case 5
        end
    end

    if case_ind==5 %case 5: all firms in top two productivity states actively selling, firms in state three randomizing, no one else actively selling
        
        %state three values reflect randomization
        th_3=th_rand;shat_3=shat_rand;S_3=S_rand;
        
        %find remaining shat's,p satisfying equilibrium conditions, given lambda 
        res_fun = @(x) resid_dim5case5_pscc(x,shat_3,S_3,la,u,mu0,phi,kap0,eta0,be,del,d,Pzz,cz,th_guess5,M);
        [x_5,fval,exit_flag]= fsolve(res_fun,[shat_guess5(1:2);al_guess5],options);exit_flag
        shat_1=x_5(1);shat_2=x_5(2);al_3=x_5(3);
        
        %find corresponding theta's
        th_res = @(t) la-mu_pscc(t,mu0,eta0)*eta_el_pscc(t,mu0,eta0)/(1-eta_el_pscc(t,mu0,eta0))*shat_1*kap_s_pscc(shat_1,kap0,phi)/(1+eta_pscc(t,mu0,eta0)*shat_1);
        [th_1,fval,exit_flag1] = fsolve(th_res,th_guess5(1),options);
        th_res = @(t) la-mu_pscc(t,mu0,eta0)*eta_el_pscc(t,mu0,eta0)/(1-eta_el_pscc(t,mu0,eta0))*shat_2*kap_s_pscc(shat_2,kap0,phi)/(1+eta_pscc(t,mu0,eta0)*shat_2);
        [th_2,fval,exit_flag2] = fsolve(th_res,th_guess5(2),options);
        
        %update initial guesses
        shat_guess5=[shat_1;shat_2];th_guess5=[th_1;th_2];al_guess5=max(al_3,0);

        if al_3>=0, case_ind=5; disp('case 5');%if randomize with positive probability in state three, label as case 5  
            %produce corresponding transition matrix and value vectors
            A=[eta_pscc(th_1,mu0,eta0)*shat_1+1,0,0,0,0;0,eta_pscc(th_2,mu0,eta0)*shat_2+1,0,0,0;0,0,al_3*eta_pscc(th_3,mu0,eta0)*shat_3+1,0,0;0,0,0,1,0;0,0,0,0,1];
            shat=[shat_1;shat_2;shat_3;0;0];th=[th_1;th_2;th_3;0;0];
            al=[1;1;al_3;0;0];%probabilities of active selling across states
            Pr=(u-be*(1-d)*(1-del)*la)/(1-be*(1-d)*(1-del))*ones(n,1);
            for ii=1:n, Ps(ii,1)=Pr(ii,1)-la/mu_pscc(th(ii,1),mu0,eta0);end;
            pr=Pr-be*(1-d)*(1-del)*Pzz*(Pr+diag(al)*(Ps-Pr));ps=Ps-be*(1-d)*(1-del)*Pzz*(Pr+diag(al)*(Ps-Pr));
        else case_ind=4;%otherwise move to consider case 4
        end
    end

    if case_ind==4 %case 4: all firms in top two productivity states actively selling, no one else actively selling
        
        %find shat's satisfying equilibrium conditions, given lambda 
        res_fun = @(x) resid_dim5case4_pscc(x,la,u,mu0,phi,kap0,eta0,be,del,d,Pzz,cz,th_guess4,M);
        [x_4,fval,exit_flag] = fsolve(res_fun,shat_guess4,options);exit_flag
        shat_1=x_4(1);shat_2=x_4(2);
        
        %find corresponding theta's
        th_res = @(t) la-mu_pscc(t,mu0,eta0)*eta_el_pscc(t,mu0,eta0)/(1-eta_el_pscc(t,mu0,eta0))*shat_1*kap_s_pscc(shat_1,kap0,phi)/(1+eta_pscc(t,mu0,eta0)*shat_1);
        [th_1,fval,exit_flag1] = fsolve(th_res,th_guess4(1),options);
        th_res = @(t) la-mu_pscc(t,mu0,eta0)*eta_el_pscc(t,mu0,eta0)/(1-eta_el_pscc(t,mu0,eta0))*shat_2*kap_s_pscc(shat_2,kap0,phi)/(1+eta_pscc(t,mu0,eta0)*shat_2);
        [th_2,fval,exit_flag2] = fsolve(th_res,th_guess4(2),options);

        %update initial guesses
        shat_guess4=[shat_1;shat_2];th_guess4=[th_1;th_2];
        
        %firm value from active selling versus not in state two
        test_2=-kap_n_pscc(shat_2,kap0,phi)-la/mu_pscc(th_2,mu0,eta0);

        if test_2>=0, case_ind=4; disp('case 4');%if prefer active selling in state two, label as case 4 
            %produce corresponding transition matrix and value vectors
            A=[eta_pscc(th_1,mu0,eta0)*shat_1+1,0,0,0,0;0,eta_pscc(th_2,mu0,eta0)*shat_2+1,0,0,0;0,0,1,0,0;0,0,0,1,0;0,0,0,0,1];
            shat=[shat_1;shat_2;0;0;0];th=[th_1;th_2;0;0;0];
            al=[1;1;0;0;0];%probabilities of active selling across states
            Pr=(u-be*(1-d)*(1-del)*la)/(1-be*(1-d)*(1-del))*ones(n,1);
            for ii=1:n, Ps(ii,1)=Pr(ii,1)-la/mu_pscc(th(ii,1),mu0,eta0);end
            pr=Pr-be*(1-d)*(1-del)*Pzz*(Pr+diag(al)*(Ps-Pr));ps=Ps-be*(1-d)*(1-del)*Pzz*(Pr+diag(al)*(Ps-Pr));
        else case_ind=3;%otherwise move to consider case 3
        end
    end

    if case_ind==3 %case 3: all firms in top productivity state actively selling, firms in state two randomizing, no one else actively selling
        
        %state two values reflect randomization
        th_2=th_rand;shat_2=shat_rand;S_2=S_rand;
        
        %find remaining shat,p satisfying equilibrium conditions, given lambda 
        res_fun = @(x) resid_dim5case3_pscc(x,shat_2,S_2,la,u,mu0,phi,kap0,eta0,be,del,d,Pzz,cz,th_guess3,M);
        [x_3,fval,exit_flag] = fsolve(res_fun,[shat_guess3;al_guess3],options);exit_flag
        shat_1=x_3(1);al_2=x_3(2);

        %find corresponding theta
        th_res = @(t) la-mu_pscc(t,mu0,eta0)*eta_el_pscc(t,mu0,eta0)/(1-eta_el_pscc(t,mu0,eta0))*shat_1*kap_s_pscc(shat_1,kap0,phi)/(1+eta_pscc(t,mu0,eta0)*shat_1);
        [th_1,fval,exit_flag1] = fsolve(th_res,th_guess3(1),options);

        %update initial guesses
        shat_guess3=shat_1;th_guess3=th_1;al_guess3=max(al_2,0);

        if al_2>=0, case_ind=3; disp('case 3');%if randomize with positive probability in state two, label as case 3
            %produce corresponding transition matrix and value vectors
            A=[eta_pscc(th_1,mu0,eta0)*shat_1+1,0,0,0,0;0,al_2*eta_pscc(th_2,mu0,eta0)*shat_2+1,0,0,0;0,0,1,0,0;0,0,0,1,0;0,0,0,0,1];
            shat=[shat_1;shat_2;0;0;0];th=[th_1;th_2;0;0;0];
            al=[1;al_2;0;0;0];%probabilities of active selling across states
            Pr=(u-be*(1-d)*(1-del)*la)/(1-be*(1-d)*(1-del))*ones(n,1);
            for ii=1:n, Ps(ii,1)=Pr(ii,1)-la/mu_pscc(th(ii,1),mu0,eta0);end
            pr=Pr-be*(1-d)*(1-del)*Pzz*(Pr+diag(al)*(Ps-Pr));ps=Ps-be*(1-d)*(1-del)*Pzz*(Pr+diag(al)*(Ps-Pr));
        else case_ind=2;%otherwise move to consider case 2
        end
    end

    if case_ind==2 %case 2: all firms in top productivity state actively selling, no one else actively selling
        
        %find shat satisfying equilibrium conditions, given lambda 
        res_fun = @(x) resid_dim5case2_pscc(x,la,u,mu0,phi,kap0,eta0,be,del,d,Pzz,cz,th_guess2,M);
        [x_2,fval,exit_flag] = fsolve(res_fun,shat_guess2,options);exit_flag
        shat_1=x_2(1);

        %find corresponding theta
        th_res = @(th) la-mu_pscc(th,mu0,eta0)*eta_el_pscc(th,mu0,eta0)/(1-eta_el_pscc(th,mu0,eta0))*shat_1*kap_s_pscc(shat_1,kap0,phi)/(1+eta_pscc(th,mu0,eta0)*shat_1);
        [th_1,fval,exit_flag1] = fsolve(th_res,th_guess2(1),options);

        %update initial guesses
        shat_guess2=shat_1;th_guess2=th(1);
        
        %firm value from active selling versus not in state one      
        test_1=-kap_n_pscc(shat_1,kap0,phi)-la/mu_pscc(th_1,mu0,eta0);
     
        if test_1>=0, case_ind=2; disp('case 2');%if prefer active selling in state one, label as case 2    
            %produce corresponding transition matrix and value vectors
            A=[eta_pscc(th_1,mu0,eta0)*shat_1+1,0,0,0,0;0,1,0,0,0;0,0,1,0,0;0,0,0,1,0;0,0,0,0,1];
            shat=[shat_1;0;0;0;0];th=[th_1;0;0;0;0];
            al=[1;0;0;0;0];%probabilities of active selling across states
            Pr=(u-be*(1-d)*(1-del)*la)/(1-be*(1-d)*(1-del))*ones(n,1);
            for ii=1:n, Ps(ii,1)=Pr(ii,1)-la/mu_pscc(th(ii,1),mu0,eta0);end
            pr=Pr-be*(1-d)*(1-del)*Pzz*(Pr+diag(al)*(Ps-Pr));ps=Ps-be*(1-d)*(1-del)*Pzz*(Pr+diag(al)*(Ps-Pr));
        else case_ind=1;%otherwise move to consider case 2
        end
    end

    if case_ind==1 %case 1: all firms in top productivity state randomizing, no one else actively selling
        
        %state one values reflect randomization
        th_1=th_rand;shat_1=shat_rand;S_1=S_rand;

        %find p satisfying equilibrium conditions, given lambda 
        res_fun = @(x) resid_dim5case1_pscc(x,shat_1,S_1,la,u,mu0,phi,kap0,eta0,be,del,d,Pzz,cz,M);
        [x_1,fval,exit_flag] = fsolve(res_fun,al_guess1,options);
        al_1=x_1;

        %update initial guess
        al_guess1=max(al_1,0);

        if al_1>=0, case_ind=1; disp('case 1');%if randomize with positive probability in state one, label as case 1
            %produce corresponding transition matrix and value vectors
            A=[al_1*eta_pscc(th_1,mu0,eta0)*shat_1+1,0,0,0,0;0,1,0,0,0;0,0,1,0,0;0,0,0,1,0;0,0,0,0,1];
            shat=[shat_1;0;0;0;0];th=[th_1;0;0;0;0];
            al=[al_1;0;0;0;0];%probabilities of active selling across states
            Pr=(u-be*(1-d)*(1-del)*la)/(1-be*(1-d)*(1-del))*ones(n,1);
            for ii=1:n, Ps(ii,1)=Pr(ii,1)-la/mu_pscc(th(ii,1),mu0,eta0);end
            pr=Pr-be*(1-d)*(1-del)*Pzz*(Pr+diag(al)*(Ps-Pr));ps=Ps-be*(1-d)*(1-del)*Pzz*(Pr+diag(al)*(Ps-Pr));
        else disp('No solution recovered');%no solution recovered
        end
    end
end

%check remaining equilibrium condition requiring stationarity
Mlom=Pzz'*((1-d)*(1-del)*A+d*nu*eye(5));%law of motion for customer base by state
det_val=det(Mlom-eye(5));%if determinant zero, stationary distribution exists for this lambda (if not, need to adjust lambda) 

%find stationary distribution
[evec,eval]=eig(Mlom);
Ntilde=evec(:,2);
a=sum(Ntilde)+sum(al.*th.*shat.*Ntilde);
N=Ntilde/a;
G_ad=(1-del)*(ones(n,1)+eta_pscc(th,mu0,eta0).*shat);%growth of actively selling firms
G_nad=(1-del);%growth of rest of firms
disp('growth rate among continuing firms (annual)') 
(((al.*G_ad+(1-al).*G_nad)'*N)/sum(N))^12
disp('probability of a sale')
al'*N/sum(N)
disp('average sale discount')
((pr-ps)./pr)'*(al.*N)/sum(al.*N)

%simulation for figure 1
tmax=1000;t_vec=1:tmax;n_vec(1,1)=1;
rng(1);rnd_z=rand(1,tmax+1);rnd_al=rand(1,tmax+1);%draw shocks
ind_vec(1)=1;%start from low state
for i=1:tmax
    if rnd_al(i)<=al(ind_vec(i)) 
        n_vec(1,i+1)=G_ad(ind_vec(i))*n_vec(1,i);p_vec(1,i)=ps(ind_vec(i));shat_vec(1,i)=shat(ind_vec(i));
    else n_vec(1,i+1)=G_nad*n_vec(1,i);p_vec(1,i)=pr(ind_vec(i));shat_vec(1,i)=0;
    end
    if rnd_z(1,i+1)<Pzz(ind_vec(i),1), ind_vec(i+1)=1;
    elseif rnd_z(1,i+1)<Pzz(ind_vec(i),1)+Pzz(ind_vec(i),2), ind_vec(i+1)=2;
    elseif rnd_z(1,i+1)<Pzz(ind_vec(i),1)+Pzz(ind_vec(i),2)+Pzz(ind_vec(i),3), ind_vec(i+1)=3;
    elseif rnd_z(1,i+1)<Pzz(ind_vec(i),1)+Pzz(ind_vec(i),2)+Pzz(ind_vec(i),3)+Pzz(ind_vec(i),4), ind_vec(i+1)=4;
    else ind_vec(i+1)=5;
    end
end

%select section for figure and plot
ind_init=31;
n_vec_fig=n_vec(12*ind_init:12*ind_init+7*12)/n_vec(12*ind_init);
p_vec_fig=p_vec(12*ind_init:12*ind_init+7*12);
shat_vec_fig=shat_vec(12*ind_init:12*ind_init+7*12);
t_vec_fig=1:length(n_vec_fig);

figure; subplot(1,3,1);set(gca,'FontSize',14);
plot(t_vec_fig/12,log(n_vec_fig),'k','LineWidth',1.5);xlabel('t (years)');title('A. Customer base');xlim([0 7]);xticks([0 1 2 3 4 5 6 7]);
subplot(1,3,2);set(gca,'FontSize',14);plot(t_vec_fig/12,p_vec_fig,'k','LineWidth',1.5);xlabel('t (years)');title('B. Prices: regular and sale');xlim([0 7]);xticks([0 1 2 3 4 5 6 7]);hold on; ylim([0.7 1.05]);
subplot(1,3,3);set(gca,'FontSize',14);plot(t_vec_fig/12,shat_vec_fig,'k','LineWidth',1.5);xlabel('t (years)');title('C. Selling intensity');xlim([0 7]);xticks([0 1 2 3 4 5 6 7]);ylim([0 .8]); 
set(gcf,'PaperPositionMode','auto');
%print -deps 'C:\Users\leena\Dropbox\pricing\PricingwCC\JPEmacro\publication_files\figure1.eps'

%produce figure 2: teaser pricing

%find steady state
eta_shat=(1+g)/(1-del)-1;
shat_fun=@(t) eta_shat/eta_pscc(t,mu0,eta0);
res_fun=@(t) (1-eta_el_pscc(t,mu0,eta0))/eta_el_pscc(t,mu0,eta0)*(u-c-be*(1-d)*(1-del)*kap_pscc(shat_fun(t),kap0,phi)-(1-be+be*ent)*(1-d)*(1-del)*shat_fun(t)*kap_s_pscc(shat_fun(t),kap0,phi)/(1-ent-(1-d)*(1-del)))/(shat_fun(t)*kap_s_pscc(shat_fun(t),kap0,phi)*(1-be*(1-d)*(1-del)*(1-mu_pscc(t,mu0,eta0))))*(1-ent-(1-d)*(1-del))/((1-d)*(1-del))-1;
th_tea_ss=fsolve(res_fun,.5);
shat_tea_ss=shat_fun(th_tea_ss);

%set value of lambda, confirm below that coincides with stationary equilibrium
la=0.0518134;

%initial guesses use steady state values
shat_guess_tea=[shat_tea_ss;shat_tea_ss;shat_tea_ss;shat_tea_ss;shat_tea_ss];

%continuation in standard deviation of shock
vol_vec=.001:.001:.1;
for i=1:length(vol_vec),i/length(vol_vec)
    %productivity vector
    sigz=vol_vec(i);
    phiz=sigz*sqrt(n-1);
    z=ones(n,1)+[phiz;phiz/2;0;-phiz/2;-phiz];
    cz=c*flip(z);%cost vector

    %find shat's satisfying equilibrium conditions, given lambda 
    res_fun = @(xx) resid_dim5teaser_pscc(xx,la,u,mu0,phi,kap0,eta0,be,del,d,Pzz,cz);
    [x_t,fval,exit_flag] = fsolve(res_fun,shat_guess_tea,options);exit_flag 
    shat_1_tea=x_t(1);shat_2_tea=x_t(2);shat_3_tea=x_t(3);shat_4_tea=x_t(4);shat_5_tea=x_t(5);
    
    %find corresponding theta's
    th_1_tea=(kap_s_pscc(shat_1_tea,kap0,phi)/la)^(1/(1+eta0));
    th_2_tea=(kap_s_pscc(shat_2_tea,kap0,phi)/la)^(1/(1+eta0));
    th_3_tea=(kap_s_pscc(shat_3_tea,kap0,phi)/la)^(1/(1+eta0));
    th_4_tea=(kap_s_pscc(shat_4_tea,kap0,phi)/la)^(1/(1+eta0));
    th_5_tea=(kap_s_pscc(shat_5_tea,kap0,phi)/la)^(1/(1+eta0));

    %update initial guesses
    shat_guess_tea=[shat_1_tea;shat_2_tea;shat_3_tea;shat_4_tea;shat_5_tea];

    %produce corresponding transition matrix and value vectors
    A_tea=[eta_pscc(th_1_tea,mu0,eta0)*shat_1_tea+1,0,0,0,0;0,eta_pscc(th_2_tea,mu0,eta0)*shat_2_tea+1,0,0,0;0,0,eta_pscc(th_3_tea,mu0,eta0)*shat_3_tea+1,0,0;0,0,0,eta_pscc(th_4_tea,mu0,eta0)*shat_4_tea+1,0;0,0,0,0,eta_pscc(th_5_tea,mu0,eta0)*shat_5_tea+1];
    shat_tea=[shat_1_tea;shat_2_tea;shat_3_tea;shat_4_tea;shat_5_tea];th_tea=[th_1_tea;th_2_tea;th_3_tea;th_4_tea;th_5_tea];
    p_tea=[1;1;1;1;1];
    Pe=(u-be*(1-d)*(1-del)*la)/(1-be*(1-d)*(1-del))*ones(n,1);
    for ii=1:n, Pn(ii,1)=Pe(ii,1)-la/mu_pscc(th_tea(ii,1),mu0,eta0);end
    pe=Pe-be*(1-d)*(1-del)*Pzz*Pe;pn=Pn-be*(1-d)*(1-del)*Pzz*Pe;
end

%check remaining equilibrium condition requiring stationarity
Mlom_tea=Pzz'*((1-d)*(1-del)*A_tea+d*nu*eye(5));%law of motion for customer base by state
det_val=det(Mlom-eye(5));%if determinant zero, stationary distribution exists (if not, need to change lambda) 

%simulation for figure 2
G_ad_tea=(1-del)*(1+eta_pscc(th_tea,mu0,eta0).*shat_tea);
tmax=1000;t_vec=1:tmax;n_vec(1,1)=1;
rng(1);rnd_z=rand(1,tmax+1);
ind_vec(1)=1;%start from low state
for i=1:length(t_vec)
    n_vec(1,i+1)=G_ad_tea(ind_vec(i))*n_vec(1,i);pe_vec(1,i)=pe(ind_vec(i));pn_vec(1,i)=pn(ind_vec(i));shat_vec(1,i)=shat_tea(ind_vec(i));
    if rnd_z(1,i+1)<Pzz(ind_vec(i),1), ind_vec(1,i+1)=1;
    elseif rnd_z(1,i+1)<Pzz(ind_vec(i),1)+Pzz(ind_vec(i),2), ind_vec(1,i+1)=2;
    elseif rnd_z(1,i+1)<Pzz(ind_vec(i),1)+Pzz(ind_vec(i),2)+Pzz(ind_vec(i),3), ind_vec(1,i+1)=3;
    elseif rnd_z(1,i+1)<Pzz(ind_vec(i),1)+Pzz(ind_vec(i),2)+Pzz(ind_vec(i),3)+Pzz(ind_vec(i),4), ind_vec(1,i+1)=4;
    else ind_vec(1,i+1)=5;
    end
end

%select section for figure and plot
ind_init=31;
n_vec_fig=n_vec(12*ind_init:12*ind_init+7*12)/n_vec(12*ind_init);
p_vec_fig=p_vec(12*ind_init:12*ind_init+7*12);
pe_vec_fig=pe_vec(12*ind_init:12*ind_init+7*12);pn_vec_fig=pn_vec(12*ind_init:12*ind_init+7*12);
shat_vec_fig=shat_vec(12*ind_init:12*ind_init+7*12);
t_vec_fig=1:length(n_vec_fig);

figure; subplot(1,3,1);set(gca,'FontSize',14);
plot(t_vec_fig/12,log(n_vec_fig),'k','LineWidth',1.5);xlabel('t (years)');title('A. Customer base');xlim([0 7]);xticks([0 1 2 3 4 5 6 7]);ylim([0 .205]);  
subplot(1,3,2);set(gca,'FontSize',14);plot(t_vec_fig/12,pe_vec_fig,'k','LineWidth',1.5);xlabel('t (years)');title('B. Prices: regular and teaser');xlim([0 7]);xticks([0 1 2 3 4 5 6 7]);hold on;
plot(t_vec_fig/12,pn_vec_fig,'k','LineStyle','-.','LineWidth',1.5);
subplot(1,3,3);set(gca,'FontSize',14);plot(t_vec_fig/12,shat_vec_fig,'k','LineWidth',1.5);xlabel('t (years)');title('C. Selling intensity');xlim([0 7]);xticks([0 1 2 3 4 5 6 7]); 
set(gcf,'PaperPositionMode','auto');
%print -deps 'C:\Users\leena\Dropbox\pricing\PricingwCC\JPEmacro\publication_files\figure2.eps'
